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This paper describes matrix methods that have been developed for cal- 
culating compressible flow on a blade-to-blade surface of revolution. The 
methods have been fully tested to date only for the design of plane cas- 
cades to prescribed blade surface distributions; the methods will be illus- 
trated here for that problem only. Similar methods are presently being 
applied to both the direct and indirect problems and for flow on arbitrary 
surfaces of revolution in annular cascades with stream sheet thickness 
variations. It is believed that by such methods, both the direct and in- 
direct calculations can be reduced to about 60 to 90 seconds of computing. 


The trend in compressor and turbine design is toward fewer and more 
highly loaded stages. To do this and maintain high efficiency demands the 
ability to calculate in ever-increasing detail the gas flow through such a 
machine. So complex are the equations governing the flow and the 
geometries involved that practicable solutions can be found only after 
making simplifying assumptions. The degree of approximation is always a 
compromise between a realistic description of the physical processes and 
a mathematical model that can be solved within reasonable time and cost. 
This has led to design procedures which treat the flow in two stages — a 
two-dimensional through-flow calculation which neglects circumferential 
variations, followed by a two-dimensional blade-to-blade calculation in 
which the flow is assumed to take place on a surface of revolution. Al- 
though fully three-dimensional calculations are being attempted, these 
are slow and costly and have a long w r ay to go before they become design 
tools. 

It seems likely, therefore, that for a few years to come, two-dimensional 
approaches will remain the basis of most design work and, for this reason, 
it is worthwhile to make these calculations as realistic and fast as possible. 
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Methods will be described here for calculations on blade-to-blade 
surfaces of revolution. These methods are being applied both to the direct 
problem of calculating blade surface velocities when the blade geometries 
are prescribed and to the indirect problem of calculating the blade 
geometry when the blade surface velocities are prescribed. The methods 
will be illustrated by discussing the indirect problem for compressible 
flow in a plane cascade. This has been chosen because it is the only problem 
for which the methods have been fully tested to date and because the 
authors have seen no other fully compressible solution to this problem. 
It is believed that the methods described here extend easily to both the 
direct and indirect problems on surfaces of revolution with stream sheet 
thickness variations. 


MATHEMATICAL ANALYSIS 


Assumptions 

The following assumptions have been made. 

(1) The flow is steady, inviscid and irrotational. 

(2) The fluid is a perfect gas. 

(3) The total temperature is uniform across the entry to the cascade. 

(4) The flow is plane two-dimensional flow and the normal com- 
ponent of velocity is zero on the blade surface. 

(5) The cascade contains an infinite number of equally spaced blades 
of infinite length. 

The assumption of irrotationality, together with the finite difference 
approximations to the differential equations and the boundary-value 
approach to the solution of the finite difference equations, tacitly assume 
that the flow is everywhere subsonic. However, the method will formally 
produce answers with supersonic patches and, where these are small and 
the peak Mach numbers only a little above sonic, these solutions are 
probably realistic. 


Equations of Motion 

In the analysis that follows, x and y are Cartesian coordinates with x 


measured in the 


recuoii aiiu y m the 


nV>TTT1 ^ 

VV U3L. 


direct 


as 


shown in figure 1. Velocities and density are normalized with respect to 
the stagnation sound speed and stagnation density, respectively. 

The equations governing the flow are those of irrotationality and 
continuity which are, respectively 


dV x dV y 
dy dx 


( 1 ) 
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f (p^)+f (py„)= o (2) 

dx dy 

Density is related to velocity through Bernoulli’s equation 

( y-1 V /(Y - 1> 

P=|l-' Y~ {V X ‘-tVy‘) I 

Equations (1) and (2) may be satisfied identically by a potential 
function <f> and stream function defined by 





d4 ' 
dx 


= ~PV V 


It will be convenient also to work in terms of the net velocity V and 
flow direction 6, related to V x and V v by the equations 

V X =V cose 

V v = V sin d 


F 



Figure 1. — One strip of the cascade in the physical plane. 
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If we now use <f> and \f/ as independent variables instead of x and y, 
equations (1) and (2) become 


dV dd 

pV — - — F 2 — = 0 
d\p dtp 


( 3 ) 


(pF) + (pF) 2 ~ = ° (4) 

dtp dip 

and Bernoulli’s equation is 

r 7 _i -|i/(r-D 

At this stage, Stanitz (ref. 1) linearized equations (3) and (4) by ap- 
proximating equation (5) by 

1 

p ~ vr+F 2 

At the equivalent stage in the direct problem, other workers have 
arranged the equations either in the form of a pseudo Poisson’s equation, 
collecting the terms describing incompressible effects on the left in the 
form of a Laplacian and the terms describing compressible effects on the 
right in the form of a source term; or they have arranged the equations 
in the form of a general partial differential equation in which the coeffi- 
cients contained derivatives of the density p. Finite difference and 
singularity methods have then been used to solve the equations in these 
forms iteratively by guessing the source term or coefficients, solving as 
though the equations were linear, and then re-estimating the terms that 
had been guessed. Iterative methods based on these forms of arrangements 
of the equations converge slowly at high Mach numbers because the 
guessed terms are In’ no means small perturbations and important con- 
tributions are left “trailing” one cycle behind in the iterations. 

In order to introduce compressibility effects quickly into an iterative 
method, the authors consider it better to use Bernoulli’s equation to 
express the derivatives of p in terms of those of the dependent variable 
and then to collect together all terms containing any particular derivative 
of that variable. The coefficients of these variables then do not contain 
derivatives of p which have to be guessed. For the indirect cascade problem 
considered here, the term d(pV)/d<p in equation (4) should not be ex- 
pressed as p(dF/d<£) + V(dp/d<p) with p and dp/d<t> being guessed. Instead, 
equation (5) should be used to obtain 


d{ P V)=p 


/ l — [(?+ 1)/2]F 2 N 
\1-[(T— 1)/2]F 2 / 
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so that equations (3) and (4) become 


p dV dd 
V dif> d<t> 


( 6 ) 


J_ / i — \ ^ \— = c\ 

pF\i-[( 7 -i)/2]FV d<t> 


r i \ /olT79\ 


(7) 


If one second-order equation was to be obtained by eliminating be- 
tween (6) and (7), then, again, the derivatives of p introduced should be 
expressed in terms of those of F. For this problem there is, however, a 
neater approach. Define F and H by 



( 8 ) 


dU _ 1 / ! — [(7+ 1 )/2]F« \ 
pF \1 — [( 7 — 1)/2]F 2 / 


dV 


so that equations (6) and (7) become 


d JL+ s ±-o 

d<f> d\ (/ 


(9) 


(10) 


d\f/ d<j> 


(ID 


Using equation (5), equation (8) may be integrated directly for some 
values of 7. Taking 7 = f and writing z= F 2 /6, we have 

3 z 3 z 2 z 3 

F(F)=logF--+— - (12) 

Taking 7 = x an d writing z 2 = 1 — ( F 2 /5) , we have 

z 5 z 3 /l+z\ 23 

F(F) = log F+- + -+Z- log (13) 


In each case, the constant of integration has been chosen such that 
F(V)— * log F as F— >0. 

The function F will now be taken as the dependent variable and 
equation (10) written in the form 


dH af dd 
dF d<(> d\p 


(14) 


where, from (8) and (9) , 
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dH_ 1 / l-[(7+l)/2]F 2 \ 
dF p 2 Vi-[(7-1)/2]FV 

The Potential and Deflection Conditions 

In figure 1, which shows one strip of the cascade, AB, CE, GF and KH 
are dividing streamlines. At any point ( x,y ) on AB, the flow conditions 
are the same as at the point (x,y+t) on JH and similarly for CD and 
GF. The lines KA and FE are far upstream and downstream of the 
cascade, where flow conditions are uniform. Figure 2 shows the same 
diagram mapped into the -plane with A, B, C, D, and E chosen as 
\p = 0. From the definitions of 4> and \p, it follows that 

d<p=V ds (15a) 

d\p = pV dn (15b) 

Define 

A <j>L = <t>H~<t>B 
A<t>T = <t>G — <f>C 

A — 

A <f>p = 4 >g~ 4>h 

Clearly 

A <t>L = <t>J — <t>K 
A <pT = <t>B — <t>D 

From (15a), remembering that J, K, E and D are far from the cascade, 





Figure 2.— One strip of the cascade in the {4>,F) plane. 
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A <j>L — V u t sin 8 U 
A<t> T = Vd t sin 6 d 

It is also clear that 

A<j>s — A if>p = Ai pL — A</>r 

so that 

Alps— A<t> P = t(V u sin 6 U —V d sin 6 d ) (16) 

Equation (16) will be called the potential condition. 

From equation (11), we have 


® F dtp+d d\p = 0 

J AEFK 

from which it follows that 
f c f° 

1 F d<p~ / F dip — \po(d u —dd)-\-F u AipL — F d Atp T (17) 

J b J H 

where ipa is the value of \p along KF . Equation (17) will be called the 
deflection condition. 

In the indirect problem, the velocity on the blade surfaces is prescribed 
as a function of fractional arc length S' measured from S' = 0 at the 
leading edge stagnation point and S' = 1 at the trailing edge stagnation 
point. Let these velocity distributions be Vs{S') and Vp(S') along the 
suction and pressure surfaces. If L s and L P represent the physical lengths 
of these surfaces, measured between stagnation points, the potential and 
deflection conditions may be written 

Ls j V s dS' — Lp ( V P dS' = t(V u sin 0„— V d sin 6 d ) (18) 

•'o •'o 

and 

Ls f (VF)sdS' Lp f (VF) P dS' = ipo(d u — 8 d ) -\-F u AipL — F d Atpr (19) 

•'o •'o 

From the prescribed velocity distributions and upstream and down- 
stream conditions, the corresponding values of F may be found from 
(12) and (13) and L s and L P from (18) and (19). The lengths in the 
(0,i^) -plane, A <p s and A <p P , may then be found from 
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and the diagram of the {<p,\p) -plane constructed. Eliminating 6 between 
equations (11) and (14) gives 


a*F d_/dH o 
dp 2 d<t> \dF d<f>) 


(20) 


To determine the blade shape corresponding to the prescribed surface 
velocity distributions and far upstream and downstream conditions, we 
have to solve equation (20) inside and on the contour ADFJ, subject to 
the boundary conditions : 


(1) F is prescribed on BC, HG, AJ and DF 

(2) Along AB and JH 

F(p,0) —F(<t>+A<t>L, Po) (21a) 

6(4>,0) Pa) (21b) 

(3) Along CD and GF 

F(<)>,0) =F(<t>-\- Apr, Po) (21c) 

6(<t>,0) =d(<t>-\-A<t>T) Pa) (21d) 


Transformation of the ( <p , ^)-Plane 

There are a number of possible approaches to a numerical solution of 
this boundary-value problem. The one given here involves an approximate 
transformation of the (p,p) -plane and some tedious algebra. However, 
the error in the transformation can be controlled so that it is less than that 
involved in the numerical methods and leads to a boundary-value problem 
posed in a form for which this is a quick and elegant method of solution. 

First, in order to get a good spacing of points on a finite difference grid 
and not to map part of the suction surface twice, it is convenient to invert 
the diagram in the (<£,^)-plane through a transformation p— >pa~ P- We 
can achieve this without altering the equations if we make the additional 
transformation 0 — * — 0. In what follows, this transformation will be as- 
sumed to have been made. Define new variables 4>' and \p' through the 
equations 

(22a) 

a 

^ = ^ + A^r ^ (ai+a2tanh ^l (22b) 

2 0 1 a J 


where 
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ft Aifty-f-Afo, 
A <f>p 2 


ft A <jyp — A <I>l 

& 2 

A dip 2 

The constant a is merely a scaling factor which can be chosen freely; 
ft is a constant which, for values of 0' > ft/2, makes tanh0'~l. T his 
transformation approximately maps the contour ADFJ of the ( <j>,\p )- 
plane into three rectangular regions in the -plane as shown in 

figure 3. If we write 

a , ft/ A <f>p 

1 + i'l'/M sech 2 <j>’ 


— (ai-f-a 2 tanh </>') 
^o[l + (^/^o)a 2 sech 2 <t>'2 

then we have 



Writing H for dH/dF, equation (20) becomes 


I 

V 



Figure 3 .—One strip of the cascade in the (<t>',F) plane. 
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. d i F dF r d . • . di/ a db' 

(a ' 2ff+6 ' 2) L a a*' ( )+ 


2b' a d~F c? d 2 F 

+ — TTT-. + — — = 0 


< 54 >' d ^ &> 2 a ^' 2 


( 23 ) 


Although this equation looks more complicated than (20), the boundary 
conditions (21a) and (21b) are simplified to 

F (</>', 0)=F(*» 

0 (<*>', O)=0(*» 

along AB and JH and similarly along CD and FG. 


Numerical Analysis 

Equation (23) may now be solved numerically by finite differences on 
a rectangular grid in the (</>', ifO -plane. The method will be described for a 
grid with spacing 8<t>' and 5^' constant in the </>' and \p' directions, respec- 
tively. In practice, it is better to use an unequally spaced grid but, to 
avoid unnecessary complication in the description, a discussion of un- 
equal grid spacing will be left until later. The grid described here is shown 
in figure 3. Write equation (23) in the form: 

»' %+* *' **' w~° (24 > 

The method of solution will be to estimate the coefficients A, B, C and 
D, solve (24) as a linear equation, and re-estimate these coefficients. The 
process is continued until converged, which usually requires about three 
or four cycles. Equation (24) may be approximated by finite differences 
in the form: 

Aj* (F/+ 1 - 2F/+F/- 1 ) +R/(F/+ 1 - F,- 1 ) + C/ (Fj+J - F }~\ - F) + _\ + F):\) 

+Z)/(F y+1 *-2F/+F y _ 1 ‘) =0 (25) 

forl<»</— 1; l<j<J—l. 

The boundary conditions are (1) that Fo’ and Fj' are given on BC and 
HG, together with F,° and F/ for j = 0 . . .J, and (2) that F a '= Fj' and 
0o * — 6j { along AB and JH and along CD and GF (eqs. 21(a)-21(d)). 
The method of solving these finite difference equations is a slight modifica- 
tion of a method suggested to the authors by Stocker (ref. 2) . Rewrite 
equation (25) , grouping terms according to superscripts *+ 1 , i, and i— 1 . 
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[ - (Aj+Bj) F/ +1 + cy-F’+i] 

+ [DfFj-i*— 2 (Aj+Dj) 

+[Cy*F;il+ (Aj—Bj) F/-» - C'/F;-}] = o (26) 

Inside the rectangle BCGH, augment equation (26) with the equations 

F 0 i =Fo i 

Fj*=Fj' 

remembering that both F a ' and Fj' are known. In the rectangles ABHJ 
and CDFG, augment equation (26) with 

Fj=Fj' 

0o < = 0p 

The last relation must be expressed in terms of F. This could be done 
using equation (11), which implies that 

( d z\-( d jy 

W/o"W/i 

and approximating this relation by finite differences. This was tried, but 
it led to small but unacceptable errors. Instead, therefore, equation (11) 
was integrated along <f > ' = constant and the boundary conditions 00 s = 6/ 
inserted into the integral. The integral was then approximated by finite 
differences using Simpson’s rule, giving 

E Kj{Fj i+ 1 —F/- 1 ) =0 (27) 

y«=o \ a / 

where Kj= 1, 4, 2 ... 2, 4, 1. 

Therefore, inside the rectangles ABHJ and CDFG, equation (26) is 
augmented by (27) andF 0 < = F/. 

Defining F‘ to be the column vector (Fj, Fj, . . . , Fj'), equation (26), 
together with the augmenting equations, may be written in the form 

M < F <+1 +iV < F < +P’F’-i= Q* (28) 

where M *, N' and P* are square matrices and Q’ is a column vector which 
contains only zeros inside the rectangles ABHJ and CDFG and is of the 
form (P 0 *, 0, 0, 0, . . . 0, F J i ) inside the rectangle BCGH. To solve equa- 
tion (28), we begin by estimating Fj at every mesh point other than 
those along t = 0 and i=l where F is prescribed. From these estimates, 
the coefficients A, B, C, and D of (24) may be calculated at each point 
and hence the matrices M', N' and P' of (28) may be determined. We 
then look for a solution of (28) of the form 
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F i = i?*F m +t i (29) 

where the R' are square matrices and the t* are column vectors. To 
determine R* and t*, we substitute (29) into (28) and, after some re- 
arrangement, obtain 

F<= - (N i +P i R i - 1 )- 1 M i F i + l +(N i +P i R i - 1 )- 1 (Q i -P i t i - 1 ) (30) 

Comparing (29) and (30) , we obtain by inspection 

Ri=-(N i +P i R i ~ 1 )- 1 M i (31) 

V= - (tf‘+P < J2*- , )- 1 (P 1 t t - 1 -Q 4 ) (32) 

Equations (31) and (32) may be solved recursively for R' and t*, for 
l<i<7— 1, once R° and t° are known. These are obtained from the pre- 
scribed value of F°, for 

F° = E°F 1 +t° (33) 

If (33) is to be satisfied, whatever the value of F 1 , we must have 

R° = 0 
t° = F° 

Having determined 7?* and t', 0<i <7— 1, we can now solve for F every- 
where, using (29) and commencing from 

F i-i = R i-i F i +t i-i 

where F 7 is the prescribed boundary condition on i = 7. Having determined 
Fj' everywhere, the coefficients A,B,C, and D of (24) may be re-estimated 
and the process repeated until successive estimates of F everywhere 
converge to within some tolerance. In practical cases, two to four itera- 
tions are usually required, depending on the level of Mach number. 

There is a further point in the calculation of F which requires dis- 
cussion; namely, the treatment of the stagnation points, the points B, C, 
G, and 77 in figure 3. Near stagnation points, V — >0 and F — ► — ». If, 
when prescribing the velocities along BC and HG, zero velocities are 
prescribed at the stagnation points, then it is clear that the methods 
described so far cannot be applied. 

A simple and approximate method of overcoming this difficulty, which 
is equivalent to removing the stagnation points by cusping the blade, is as 
follows. At the start of each compressibility iteration, a nonzero velocity 
is assigned to the points B and 77 and another nonzero velocity to the 
points C and G. With these values, together with the other prescribed 
boundary conditions, we can now solve for F everywhere by the methods 
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already described and this solution will satisfy all the prescribed boundary 
conditions. However, for arbitrary choices of velocity at the points B and 
H and C and G, the function F is not constant at upstream and down- 
stream infinity; that is, although dF/dip' is zero there, dF/d<t>' is not zero. 
Furthermore, for given boundary conditions, the value of dF/dtf at 
upstream infinity is primarily controlled by the velocity assigned to B 
and H, and dF /dtp at downstream infinity by the velocity assigned to C 
and G. Therefore, at the start of each iteration, as well as recalculating 
the matrices M l , N l , and P*, new estimates are made of the velocities at 
B and H and at C and G to make dF / d<j>' zero at points far upstream and 
downstream. This additional change does not seriously affect the con- 
vergence of the main iteration. 

Although this is a rather crude treatment of the stagnation points, it 
does lead to accurate answers in the following sense. When 6 is calculated 
from /' , equation (11) is integrated along a streamline starting from far 
downstream where 6 is prescribed. The closeness of agreement of the 
calculated and prescribed values of 8 far upstream is one measure of the 
accuracy of the calculation. This agreement is best (about 0.2 percent for 
100 of turning) when the adjustments described have converged. 
Methods such as those of Woods (ref. 3) for dealing with singular points 
were tried but did not appear to increase the accuracy of the calculation, 
possibly because the computing grid was coarse compared with the small 
region over which the velocity is close to zero. 

brom the converged solution for F, the blade coordinates may be cal- 
culated. This is done by first integrating equation (11) along \ft’ = a/2 to 
give 6 along the center of the blade passage and then integrating equation 
(10) away from this mean line to give 8 on the blade surface. Having 
found 8, the blade coordinates are found by integrating the equations 

d<t> chf/ . 

dx = — cos 8 — — sin 8 
V pV 


, dtj> , d\l/ 
dy = — sin 0-| — -cosS 
V pv 

The integration is performed in the (tf>',\f/') -plane and commences 
from arbitrary values of x and y in the middle of the blade passage, out 
along the line </>' = constant to the blade surface and then along the blade 
surfaces, \p' = 0 and \p' = a. This path of integration avoids the necessity of 
crossing the stagnation point region. The blade shapes obtained show the 
cusps over the first and last two points on each surface and the leading and 
trailing edges are generally rounded by eye. 
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SIZE AND SPEED OF COMPUTER PROGRAM 

The methods described have been programed on an IBM 360/65 
computer. Using 40 points of each blade surface, 50 upstream and 50 
downstream points, and 11 points across the blade passage, the program 
size is 162K bytes. For a fully converged solution, three to five cycles are 
required at an average of 24 seconds per cycle. For a fully compressible 
calculation on such a large grid, the method is therefore very fast. To 
obtain this speed of computation, an unequally spaced grid has been used, 
with the grid becoming more widely spaced far upstream and downstream. 
The only change required in the methods described is to modify the finite 
difference approximations to derivatives in the obvious way. 

Sample Calculation 

The program has been tested on a number of examples, one of which, a 
NASA blade taken from reference 4, is described here. In figure 4, the 
circles and triangles represent the measured velocity distribution while 
the full line is the velocity used in the calculation. The measured outlet 
angle was changed by about 0.7° to —67.7° because the calculation cannot 
take into account viscous effects. The true and calculated blade shapes are 
shown in figure 5, where it will be seen that the agreement is generally 
good. Agreement is worst near the leading and trailing edges. The shape of 
the leading edge depends critically on the velocity distribution and this is 



Figure 4. — Velocity distribution 
of blade of reference If. 
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impossible to measure at points sufficiently close together to give accurate 
definition. Also, one of the measured points on the pressure surface has 
been ignored, for it was found that a smooth velocity distribution through 
that point did not reproduce the correct blade shape. The velocity dis- 
tribution used in this region is merely guessed to give a reasonably good 
blade shape. 


LIST OF SYMBOLS 

F A function of velocity 

H A function of velocity 

n Distance normal to a streamline 

S Distance along a streamline 

S' Fractional length along a blade surface measured between stagnation 
points 
t Pitch 
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V Velocity, normalized with respect to the stagnation sound speed 
V x The ^-component of V 
V v The ^-component of V 

x Cartesian coordinate measured in the axial direction 

y Cartesian coordinate measured in the pitchwise direction 

7 Ratio of specific heats 

6 Flow direction measured counterclockwise from the positive x 
direction 

p Density, normalized with respect to the stagnation density 
4> Potential function 
< f / Transformed potential function 
ip Stream function 
i p' Transformed stream function 

Subscripts and Superscripts 

d Far downstream 
i Index referring to the value of <j>' 

j Index referring to the value of ip' 

L Leading edge 

T Trailing edge 

a Far upstream 
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DISCUSSION 


PAY NE (Rolls-Royce): The authors are to be congratulated on 
applying a highly efficient matrix method to the solution of the boundary- 
value problem which cascade design presents in the compressible flow 
function plane. 

The technique used to solve each iterate of equation (24) (as yet un- 
published by Professor Stocker) transmits boundary-value information 
just once to the right through the vectors t and just once to the left 
through the F vectors themselves. This elegant technique thrives on 
highly rectangular grids, such as the one established here in the (</>',$')- 
plane, although the slightly approximate transformation (22b) into this 
plane could possibly be avoided by the use of a variable skew mesh in 
the (<M0 -plane. 

The starting approximation to the coefficients in equation (28), 
although not explicitly stated, presumably results from taking H = F, and 
this assumption is, in itself, quite accurate for Mach numbers less than 
about 0.8 (ref. D-l). 

The desirability of basing the design of gas turbine blading on a pre- 
scribed distribution of surface velocity can be justified by consideration 
of the mechanical, aerodynamic, and mathematical aspects of the overall 
design problem (ref. D-2). For the past eight years, all turbine blades 
designed at the Bristol Engine Division of Rolls Royce Ltd. (previously 
Bristol Siddeley Engines) have been produced on this basis, using the 
Bristol Design 1 ransformation (ref. D-3) to generate the necessary 
cascade geometry. Until now, it has been a rather unfortunate handicap 
that, while a complete velocity distribution theory (parametric descrip- 
tion of velocity distribution, optimization under geometric constraints, 
wake models, etc.) could be established for an arbitrary density-speed 
relation (ref. D-l), the actual transformation from the {<p ,$) -plane to the 
(A>2/)-plane was only practical, on a routine basis, for a simplified form of 
the p(F) function (linearized compressible flow, or Chaplygin gas). 
Although the linearized transformation can be shown to agree closely 
with plane-flow experiments for Mach numbers up to about 0.85 (ref. 
D-l), it is to be expected that the methods of Mr. Silvester and Miss 
Fitch w ill produce a routine design transformation able to cope with near- 
sonic and, perhaps, slightly supersonic flow, as well as allowing incorpora- 
tion of blade-to-blade variations of radial aerodynamic influences, as such 
variations become better understood. 
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J. P. GOSTELOW (Cambridge University) : The authors introduce 
their promising new matrix techniques as being suitable for both the 
direct and the indirect problems of cascade flow prediction. Since it is 
well known that some considerable effort at their company (ref. D-4) 
has been invested in iterative solutions of the direct problem, using the 
Martensen method (ref. D— 5) as a basis, it would be of interest to know 
whether the iterative approach has failed and has therefore been written 
off. The difficulty with the iterative schemes is that the source distribution 
contains first derivatives of p and, therefore, second derivatives of \[/. It 
would not be surprising, therefore, if convergence difficulties were experi- 
enced at high subsonic Mach numbers where the desired result is masked 
by rounding-off errors. This question does not concern simply the direct 
problem since, as Murugcsan and Railly (ref. D— 6) have shown, the 
Martensen method can become a successful design tool in solving the 
indirect problem. 

It is interesting to observe that Silvester and Fitch deliberately re- 
arrange the equations so that density-change information is transmitted 
immediately into the flow solution. It is more conventional for the density 
calculation to lag the stream function calculation by one iteration, again 
deliberately to improve stability. This latter approach is employed in 
Smith’s excellent paper (ref. D— 7) and in most streamline curvature 
solutions to the axisymmetric problem. It was clear from Smith s presenta- 
tion that the trailing density approach is justified for cases where the 
local Mach number does not exceed 0.85, but even linearized flow models 
can cope with such examples. It w'ould be interesting to know w'hether 
Smith can retain numerical stability, with density lagging by one itera- 
tion, w hen sonic conditions are reached and exceeded on the blade surface. 

The kernel of the question is whether one ought to follow Silvester and 
Fitch in rearranging the equations w'hcn local sonic conditions are ap- 
proached. 


L. MEYERHOFF (Eastern Research Group) : I have three questions. 
The first is about trailing edge conditions. I’m curious to know 
what the author believes w r ould happen if, in our reiteration, the stagna- 
tion point of the trailing edge was set right at the trailing edge to zero 
velocity and the iteration continued with that fact reinserted in each 
iteration. The other questions are (1) What is meant by the term “cycle” 
for fully converged solution in your report? Is the word cvcie meam to 
be “iteration number”? and (2) What is the total number of mesh points 
allowed by the program at present? 


H. YEH (University of Pennsylvania) : You refer to the need for an 
estimate by the computer for the velocities near the inlet and the trailing 
edge in order to have the prerequisite velocity at plus and minus infinity. 
Now, isn’t this due to the fact that they really cannot completely describe 
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the velocities anyway because you must have prerequisite separation for 
the whole profile to fulfill your inlet and exit conditions at infinity. 
Furthermore, there is a condition for which you get a closed profile. In 
other words, the profile may not be closed; this is a so-called conditioned 
closure. Now, if you had to make use of these conditions beforehand, and 
if you had considerable freedom in adjusting the velocity distribution and 
did so, it seems to me that you would not really need a computer to make 
further adjustments. 

J. W. DZIALLAS (General Electric Co.): Here are a few questions 
which should be of general interest. 

(1) If the flow is assumed everywhere subsonic, how can the field 
contain “supersonic patches”? If there are these patches, where are they 
located? Doesn’t the authors’ selection of the function F(V) near the 
stagnation points strongly affect these supersonic patches? 

(2) Is the Kutta condition satisfied? 

(3) How close to the sonic velocity can the authors’ method go on the 
profile surface? Does the solution become unstable? 

(4) What useful information can the authors extract from their 
hodograph? 

(5) Recalling the comparison with the experimental velocity dis- 
tribution presented in a slide, I ask : How valid is this comparison since, 
through smoothing of the data, adjustments on the function F(V), and 
variable grid size it seems possible to arrive at predetermined results. 
How many trials are necessary to recover the profile? 

(6) It would be interesting to see a comparison with an exact direct- 
method airfoil computation. 

P. N. R. SHEKHAR (University of Liverpool, England) : At Liverpool 
University, we have been concerned with the problem of designing airfoils 
in two-dimensional cascades. Hence, we would like to raise the following 
points : 

(1) Equation (17) is valid only for special cases of y. According to 
reference D-8, the deflection condition for any y is (fig. 2) 

I (log q d<j>+d/p <fy) — f f 0[d(l/p)/d<j>]d</># = O 

J ■ •'^*=0 oo 

This equation can only be solved iteratively. The final solution is con- 
sistent with the Price- Martensen theory (ref. D-9). 

(2) In the main paper, the problem is considered as well posed in the 
(<!>', i ') -plane. How r ever, the problem can be well posed in the {<t>,4 / ) -plane 
itself by Green’s function of the second kind as demonstrated in reference 
D-8. Hence, one wonders if it is not advantageous to work in the (<f>,^)- 
plane itself? 
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(3) Once the boundary conditions are formulated, the problem can 
be treated as a typical boundary-value problem. The methods available 
include (i) Green’s function, (ii) finite difference scheme, and (iii) varia- 
tional finite element. Even Stocker’s method (ref. 2) could be used ad- 
vantageously. However, in reference D-8, it is used to get the following 
matrix: 

[^ ][?] = [/] 

where A is a codiagonal block matrix with submatrices that are also 
codiagonal. No equation contains more than five nonzero elements and 
only the nonzero elements are stored; hence, the storage requirement is 
minimized. This method is attractive compared to the marching procedure 
for two reasons, at least. 

(a) The distribution of log q is found at all the interior points of 
the rectangle in one go. 

(b) The boundary conditions are consistent with the interior 
solution, whereas in marching procedures this is not so. In our opinion, 
once the boundary conditions are known, it really does not matter which 
method is adopted for determining log q inside the rectangle. We are sure 
Stocker’s method could also be adopted very effectively. 

(4) At stagnation points, logg has logarithmic infinity and 0 is 
multivalued. According to L. C. Woods, the movement of the front 
stagnation points by TT >Vt> of chord distance affects the velocity peak by 
more than 10 percent for isolated airfoils (let alone cascades). What is 
really important is not so much the presence of the stagnation points as 
the effect it might have on the rest of the solution. It is probably true that 
Woods’ method needs a very refined mesh. However, Payne has proposed 
a very attractive method for determining the effect of stagnation points 
on the rest of the solution by integral equation techniques. A detailed 
analysis is available in references D-8 and D-l. We find it very difficult 
to accept the concept that the solution achieved by ignoring four stagna- 
tion points is satisfactory. One could even say that the classic channel 
model proposed by Stanitz is satisfactory for cascades. Stanitz has pro- 
duced some very realistic profiles in NACA 1116. 

(5) Last, we would like to examine the following two problems: 

(a) Inconsistency with the Price-Martensen theory 

(b) The simplicity of Green’s function solution. 

The Price-Martensen theory has been used extensively and, to a great 
extent, satisfactorily. However, we understand from Silvester-Fitch that 
Smith’s solution is consistent with the design problem. Presumably this 
means that the stagnation points, shape of the stagnation, and other 
streamlines tie up completely. Hence, it should be pointed out that with 
the help of design and Smith problem, and treated on an iterative basis, 
it should be possible to produce a one-to-one correspondence and a closed 
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profile. The simplicity of Green’s function solution can be illustrated very 
simply by taking the incompressible or linearized flows. The profile shape 
depends completely on the boundary conditions and it is immaterial what 
is happening inside the boundary. In the method reported, it is necessary 
to know what is happening not only on the boundary but inside the 
boundary as well. The difference between the incompressible and com- 
pressible flow lies in the presence of a double integral term and some 
minor points. 

SILVESTER and FITCH (authors) : The authors agree with Dr. 
Payne about the desirability of being able to design blades to prescribed 
surface velocity distributions. We believe that both the direct and in- 
direct approaches can be useful to the blade designer and it was for this 
reason alone that we developed an indirect method alongside our existing 
direct method based on the work of Martensen and Price. 

We have found that the Martensen-Price method converges well for 
subsonic flow and that convergence can be obtained, although somewhat 
more slowly, for flows containing supersonic patches, provided these are 
not too large. We have also been looking at matrix methods for the direct 
problem because we believe that they can be made faster than singularity 
methods. We also believe, but have not shown, that the immediate 
introduction of density change terms into the equations will improve the 
stability and rate of convergence. 

Concerning details of the calculation, we agree with Dr. Payne that 
we could have worked on a skew mesh in the ( 0 , 1 /') -plane. With the 
method adopted, the algebra is more tedious, but this is compensated for 
by the fact that the approximation of partial derivatives by finite differ- 
ences with small truncation error and using only the lines, i— 1, i, i+1 is 
easier in the (0',0') -plane. 

We do not commence the calculation by assuming H = F. Referring to 
figure 3, w y e assume velocities of V u on AB and JH, Va on CD and GF, 
and the prescribed velocities on BC and 11 G. The velocity elsewhere is 
assumed to vary linearly with \p' at constant 0'. 

In response to the comments of Meyerhoff and Yeh, it is certainly true 
that for given values of V and 9, far upstream and downstream, not every 
velocity distribution that the designer may prescribe wall give a closed, 
nonintersecting curve for his blade profile, but only those velocity dis- 
tributions which satisfy the so-called closure conditions. It is also true 
that although we have tacitly assumed a closed profile when deriving the 
equations, we have not placed any restrictions on the velocity distributions 
that may be prescribed and so may not obtain sensible blade shapes for 
every velocity distribution. When we use this program, we assume that 
the designer is able to specify a velocity distribution which nearly satisfies 
the closure condition and which will require only slight modification within 
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the program. In practice, we use the program in conjunction with a 
direct method, so this is usually true. We justify this approach with two 
arguments. First, when designing a cooled turbine blade, there are factors 
other than aerodynamic (stressing and cooling considerations) which 
place restrictions on an acceptable blade geometry. In order to satisfy all 
these conditions, it is likely that the designer will require three or more 
runs of the program to achieve a satisfactory blade, modifying his velocity 
distribution with each successive run. Since the program prints out the 
modifications to the velocity distributions that it makes internally and 
these can be fed into the next run, after the first one or two runs, little if 
any internal modification is required. Second, we know of no method of 
determining velocity distributions for compressible flow satisfying the 
closure condition which would involve the designer in any less work than 
the method we use. 

The choice of nonzero velocities at the stagnation points is a necessity 
since it is impossible to evaluate F or H for zero velocity. 

If arbitrary, nonzero values are used and reinserted each cycle, a con- 
verged solution may or may not be found. If convergence is obtained, 
then it will be found that, although the velocity attains the values V u 
and V d far upstream and downstream, dV/d<t> is not zero there. In addi- 
tion, integration of dQ/dcj) between far upstream and far downstream will 
not result in the prescribed turning, d d —d u . The authors regard the finding 
of velocities at the stagnation points as a process of finding the shape and 
direction of cusps which must be added to the rounded profile in order to 
support the prescribed velocity on the remainder of the profile. 

In the paper, the words “cycle” and “iteration” have both been used to 
denote the process of solving numerically equation (24) for fixed estimates 
of the coefficients A , B, C, and D at every point and of the velocities of 
the stagnation points. 

The program allows up to 50 points upstream and downstream of the 
blade, 40 points on each surface of the blade, and 11 points across the 
blade passage. Best results have been obtained by using a mesh of variable 
spacing in both the </>' and <// directions, having points closer together near 
the blade surfaces and particularly so near the leading and trailing edges. 

As to the remarks of Dziallas, it must be made clear that the indirect 
method described here is intended for use as a design tool in which blade 
shapes are determined from velocity distributions prescribed by the 
designer. Although we have tested the program by using it as a direct 
method (that is, by trying to recover blade shapes from measured velocity 
distributions) this is not the mode in which it was intended the program 
should be used. The method described here is not intended as an alterna- 
tive to the direct methods but as an additional aid to the blade designer. 
With this intended use of the program in mind, such questions as, “How 
many trials are necessary to recover a profile?” are not strictly applicable. 
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Taking the last two questions first, the program has been tested by 
recovering blade shapes both from measured velocity distributions and 
from distributions calculated by the singularity method of Martcnsen 
and Price, the matrix method of D. J. Smith, and the author’s own matrix 
method, which is still under development. (The method of Katsanis was 
not available to the authors.) It was found that, for a given blade, the 
three direct methods gave velocity distributions which differed by small 
but significant amounts, so that each produced a slightly different blade 
with the indirect program. Best agreement was with the authors’ own 
direct method. Because of the inconclusive nature of these tests and 
because it is a more valid test of a program to a designer, most testing of 
the program has been with measured velocity distributions. In those cases 
where the tests were two-dimensional and shock-free, blade shapes could 
be recovered well with velocity distributions close to those measured. 
There is some freedom in choosing F near the leading and trailing edges 
because pressure tappings are rarely close enough to give an adequate 
picture there. Such a test of the program is a useful one, provided one 
remains close to the measured results (which may include small experi- 
mental errors) . It simply is not true that one can arrive at predetermined 
results. Blade shapes are independent of the mesh size, provided it is 
fine enough. 

Although the numerical methods used can be justified only for elliptic 
equations and hence subsonic flows, even so we can and sometimes do 
prescribe velocity distributions with supersonic regions. If converged 
answers are obtained, these must necessarily contain supersonic patches 
adjacent to the regions of the blade when supersonic velocities have been 
specified. It can be assumed that in these supersonic regions, small errors 
due to round-off increase the more distant a mesh point is from the 
boundary where the velocity is specified. If the patches are small, ihe 
errors may not have a chance to grow too large, so it seems possible that 
sensible answers may be obtained. There is, however, no provision in the 
program for discontinuous solutions as would be caused by shocks. 

The program effectively selects F near stagnation points. We still have 
some further work to do on this, but we have found that it is more difficult 
to converge on velocities at the stagnation points when supersonic patches 
have been prescribed. 

As for the Kutta condition, remember that we produce cusped blades. 
The velocities on the cusps (that is, the velocities which are chosen to 
satisfy upstream and downstream boundary conditions) are chosen to be 
equal on both pressure and suction surfaces. In addition, we usually 
prescribe velocity distributions which become equal on both surfaces 
close to the leading and trailing edges. This treatment is something like a 
Kutta condition, although we do not talk specifically in terms of zero 
velocity. 
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We have not looked seriously at solutions in the hodograph plane. 

In reply to Mr. Shekhar, if the function F is defined by equation (8) , 
then equation (17) is true for all values of 7 and for all Mach numbers. 
Closed analytic forms for F can be found for rational values of 7, and 
the values of \ and | given by the authors should be quite sufficient for 
all practical purposes. The error involved in using 7 = 4/3 for the design 
of a hot turbine blade will be negligible compared with the errors intro- 
duced by other assumptions in the mathematical model such as, for 
example, isentropic inviscid flow. If this is accepted, then it is certainly 
more convenient and almost certainly more accurate to solve the potential 
and deflection conditions explicitly by simple numerical integration of the 
prescribed velocity distribution than iteratively by methods requiring 
the numerical evaluation of double integrals over the whole flow field 
with each compressibility cycle. Also, it is worth mentioning that the 
deflection condition, as we have formulated it, depends only on the 
boundary conditions. This is useful in two ways. First, it enables the 
surface lengths to be calculated without any knowledge of the flow field 
elsewhere and so allows the possibility of abandoning the program before 
any major computation has been performed if the prescribed velocity 
leads to unrealistic surface lengths. Second, it allows a check of the ac- 
curacy of the solution of equation (24) for fixed boundary conditions, 
because when the coordinates (x,y) of the blade are eventually found from 
the solution of (24), the lengths can be calculated and compared against 
those calculated from the potential and deflection conditions. We have 
found that for fully converged answers, the lengths calculated by the two 
methods agree to within less than 0.2 percent. 

The relative advantages of the (<f>,^) and (<)>', \ p') planes have been given 
in the reply to Dr. Payne. 

With reference to question (3) , it seems as though Mr. Shekhar believes 
that the method used by the authors is a marching procedure and one in 
which, in some way, the solution is inconsistent with the boundary condi- 
tions. We do not use a marching procedure; the solution obtained depends 
at every point upon all the boundary conditions and is completely con- 
sistent with them. Moreover, the solution is obtained “all in one go’’ just 
as much as in his own method, for (using Mr. Shekhar’s own notation), 
Stocker’s method is simply a method of solving the matrix equation 

Concerning question (4), again, it is not true to say that the authors 
have neglected the stagnation points. There is a striking similarity be- 
tween the method used by the authors and the treatment described by 
Payne as a relaxed treatment. Referring to figure 3, Payne’s method 
consists of the following steps : 
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(1) Choose nonzero velocity at the stagnation points. 

(2) Set the velocity at V u along AKi and JK 2 and V d on K^F and 
KzD. 

(3) Apply the cyclic or repeat conditions to the segments K^H, KiB, 
GK t and CK S . 

(4) Set up a variable, but one-parameter, mesh spacing. 

(5) Solve as though there were no singularities. 

Payne points out that if a solution is now obtained, ignoring the singu- 
larities and with an arbitrary mesh parameter, the points B and H, for 
example, will not be one blade pitch apart in the physical plane — but for 
a particular choice of the mesh parameter, this can be achieved. Notice 
that by making the assumptions (2), Payne is, in fact, forcing dV/d<t> = 0 
far upstream and downstream but, at the same time, relaxing the repeat 
conditions on 0 over the segments on which V is prescribed. Relaxing 
these conditions permits the streamlines JH and AB to be of different 
shape when they should be identical, but the error is reduced by forcing 
B and H to be nearly one pitch apart. 

This approach is very similar to that of the authors. We apply the 
repeat conditions on both V and 0 over the whole length of the dividing 
streamlines, so that the streamlines corresponding to AB and JH, for 
example, are identical in shape. It then follows that because A and J are 
one pitch apart in the physical plane, B and H must be. If we were to 
follow Payne and choose fixed velocities at the stagnation points, then we 
would have to choose a mesh spacing to make d V/ d <}> = 0 far upstream 
and downstream. Instead, we keep the mesh spacing fixed and vary the 
velocity at the stagnation points. The mesh spacing and the chosen 
velocity are to some extent interchangeable, for both affect the calculated 
values of derivatives of V near the leading and trailing edges. Payne also 
justifies the use of such approximate methods of dealing with stagnation 
points. 

As already stated in reply to Mr. Dzillias, the authors do not obtain 
complete agreement with any direct method, just as none of the direct 
methods is in complete agreement with any other. The differences are 
small, but some further work is needed. 

Finally, I would agree that where the equations of motion can be 
reduced to Laplace’s equation (that is, for two-dimensional incompressible 
flows or flows of sufficiently low Mach number that one could reasonably 
assume H = F) an integral method is probably to be preferred to a differ- 
ential equation approach. Most practical cases, however, cannot be 
described adequately by Laplace’s equation, either because the Mach 
number level is too high or because it is necessary to take into account 
effects such as stream-tube thickness variation or the fact that a turbine 
blade row does not form a linear two-dimensional cascade. (These effects 
have yet to be incorporated in the authors’ program.) Therefore, in most 
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practical cases, it is necessary to compute the fluid velocity everywhere 
(not only on the boundaries) , whichever method is used. In such cases, 
it is debatable whether integral methods are to be preferred to differ- 
ential methods. It should be pointed out that the double integral to which 
Mr. Shekhar refers is not simply a minor term in integral methods; apart 
from increasing the amount of computation to be done (compared with 
incompressible flow) it does express the difference between incompressible 
and compressible flow and this can be quite marked when Mach number 
levels are high. Further, if the other effects referred to were included, the 
double integral term would express the difference between plane in- 
compressible flow and compressible flow with stream-tube thickness 
variation and in an annular cascade. 
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